Enhanced One-Loop Corrections to WIMP Annihilation 
and their Thermal Relic Density in the Coannihilation 

Region 

Manuel Drees and Jie Gu | 

Physikalisches Institut and Bethe Center for Theoretical Physics, 
Universitat Bonn, Nussallee 12, 53115 Bonn, Germany 

Abstract 

We consider quantum corrections to co-annihilation processes of Weakly Interacting 
Massive Particles (WIMPs) due to the exchange of light bosons in the initial state ("Som- 
merfeld corrections"). We work at one-loop level, i.e. we assume that these corrections can 
be treated perturbatively. Co-annihilation is important if there is at least one additional 
new particle with mass close to the lightest WIMP, which is a Dark Matter candidate. In 
this case the exchange of a (relatively light) boson in the initial state can change the iden- 
tity of the annihilating particles. The corrections we are interested in factorize, as in the 
case of WIMP self-annihilation treated previously, but they can mix different tree-level 
amplitudes. Moreover, even small mass splittings between the external particles and those 
in the loop can change the relevant loop functions significantly. We find exact analytical 
expressions for these functions, and illustrate the effects by considering the cases of wino- 
or higgsino-like neutralinos as examples. 
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1 Introduction 



Decades after the discovery of the first hint of the existence of Dark Matter, still little is 
known about its nature. The perhaps best motivated Dark Matter candidates are "weakly 
interacting massive particles" (WIMPs), since they automatically have the roughly correct 
relic density if they are produced thermally within standard cosmology. These particles would 
have decoupled from the thermal plasma at a temperature of about 5% of their mass, when they 
were already quite non-relativistic. Their relic density scales inversely with the total (effective) 
WIMP annihilation cross section. Precise predictions for the relic density, with an error similar 
to or smaller than that of current or near-future determinations of the overall Dark Matter 
density from cosmological observations, therefore require precise calculations of the relevant 
annihilation cross sections. Right now the best observational constraint on the Dark Matter 
relic density comes from the combination of data from the WMAP satellite, type la supernovae, 
and the baryon acoustic oscillation (BAO), which gives [I] 

Q CDM h 2 = 0.1131 ±0.0034, (1) 

where h is the scaled Hubble parameter defined by H = 100/i km s -1 Mpc~ . Data from the 
PLANCK satellite should soon reduce the error to about 1.5%. Accordingly, on the theory 
side at least leading loop corrections to the annihilation cross section of WIMPs should be 
calculated to attain a percentage level precision. 

One class of potentially large loop corrections is due to long range interactions between 
the WIMPs before their annihilation, mediated by the exchange of a boson with mass well 
below the WIMP mass. Consider two WIMPs with mass m x coming to a head-on point-like 
annihilation. Each WIMP can be described by a plane wave function. In the classical limit the 
exchange of a boson <fi with coupling parameter a before the annihilation produces a potential, 
which is Coulomb-like if (f) is massless or Yukawa-like otherwise. If the Bohr radius l/(am x ) 
is smaller than the interaction range l/m^, the plane wave functions are significantly deformed 
within the potential. As a consequence the annihilation cross section is enhanced (suppressed) 
in the case of an attractive (repulsive) potential. The magnitude of the correction depends on 
the strength of the potential. The larger the Bohr energy of the potential a 2 m x /2 is compared 
to the kinetic energy m x v 2 /2 of the dark matter particle, the larger is the correction to the cross 
section. If both conditions are satisfied, i.e. am x /m^ > 1 (radius condition) and v < a (energy 
condition), the correction to the cross section is so strong that the perturbative expansion 
breaks down and one has to solve the Schrodinger equation with the potential to compute the 
deformation of the wave functions [21 13, S]. 

However, the radius condition am x /m ( j ) > 1 can raise naturalness issues in any extension 
of the Standard Model of particle physics that aims to ease the hierarchy problem. In the 
Minimal Supersymmetric extension of the SM (MSSM), for example, this condition can hardly 
be realized if is a weak gauge or Higgs boson. In the most natural case, with the mass of 
Xi < 1 TeV and the weak coupling aw ~ 1/30, am x /m ( f ) is still smaller (but not necessarily 
much smaller) than one. A one-loop calculation, which can be performed analytically, should 
then still produce a reasonably good approximation to the "exact" cross section based on a 
fully non-perturbative calculation. 

This paper builds on ref . [3] , where enhanced one-loop corrections to WIMP self-annihilation 
were treated. We extend these results to the co-annihilation region where the WIMP sector 
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includes several particles with masses close to that of the lightest WIMP. This, e.g., occurs nat- 
urally if the Dark Matter particle is part of a non-trivial multiplet of some non-abelian gauge 
group, which is broken at a scale somewhat smaller than the WIMP mass. This situation is 
considerably more complex than that treated in ref. [5J. On the one hand, co-annihilation 
between the lightest WIMP and one of its slightly heavier siblings has to be included, along 
with annihilation of the heavier states; the Sommerfeld-enhanced one-loop corrections will in 
general be different for these different initial states. Moreover, the exchange of a boson in the 
initial state can change the identity of the annihilating particles, as illustrated in Fig. [T] Even 
in the usual non-relativistic approximation, where the corrections factorize at the amplitude 
level, the one-loop correction to a given annihilation amplitude can therefore be proportional 
to a different tree-level annihilation amplitude. This means that the corrections no longer 
factorize at the cross section level, and signs - or, more generally, phases - between tree-level 
annihilation amplitudes with the same final state but different initial states become relevant. 




Figure 1: One loop correction to the annihilation of the WIMPs XiX2- They change into a new 
particle pair X3X4 after the exchange of boson 0. The four dark fermions ^j, (z = 1, . . . , 4) may 
or may not be different states, but we assume them to be close in mass. 

In this paper we obtain general analytical formulas for the one loop "Sommerfeld" correction, 
allowing all four masses before and after the boson exchange to be (slightly) different; this makes 
the use of the fitting formulae of ref. [5] unnecessary. We will see that the mass splitting is 
relevant unless \5m\m x <C rr&, where is the mass of the exchanged boson. A finite mass 
splitting can even enhance the size of the corrections, if the particles in the loop are somewhat 
heavier than the external particles; this agrees with results of ref. [6], where purely off-diagonal 
couplings to an "excited WIMP" were treated in the non-perturbative regime. As usual, we 
treat annihilation from S- and P-wave initial states separately; the resulting correction factors 
differ whenever the exchanged boson has a non-vanishing mass. 

Our results are applicable to any co-annihilation process in any model. Two concrete sample 
calculations are given in the framework of the MSSM, and the corresponding corrections to the 
relic density are shown. This complements refs. (2J [3, [8], where these corrections were computed 
in the non-perturbative regime, i.e. for wino-like neutralinos with masses well above 1 TeV. 

The remainder of the paper is organized as follows. Section [2] introduces the generalized 
formalism based on [S] and gives the model-independent analytical expressions for the correc- 
tion factors at the amplitude level. We also point out that the exchange of fermions, which can 
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occur if the co-annihilating WIMPs have different spin, does not lead to enhanced corrections. 
Section [3] discusses some properties of the solutions. In Section [| we apply these results to two 
MSSM scenarios. The last Section summarizes. 



2 Formalism 

Consider the general one-loop process depicted in Fig. [TJ involving the exchange of a relatively 
light boson between two fermions in the initial state and two possibly different fermions in 
the intermediate state: 

4> arm. ^ r , r ,~\ 

Xi + Xi — > X3 + X4 > X + Y. (2) 

excn. 

We will assume that all four fermions^ are close in mass to the lightest WIMP, which is a Dark 
Matter candidate. X and Y in the final state are standard model particles. We are interested in 
computing the annihilation cross sections during and after the decoupling of the WIMPs from 
the plasma of SM particles. Since decoupling occurs at temperature T ~ m x /20, an expansion 
of all relevant amplitudes in terms of the relative velocity v, or of the three-momentum p 
of the annihilating particles in the center-of-mass system (cms), in most cases [9] converges 
reasonably fast. The annihilation amplitude can be decomposed into partial waves, and only 
the leading S- and P-wave contributions are important, which start at order \p\° and Ipp, 
respectively. 

We are interested in scenarios where the the boson mass is significantly smaller than the 
typical WIMP mass m x . The dominant contribution to the loop amplitude then comes from 
configurations where the virtual momentum carried by <fi is much smaller than the momentum 
exchanged in the X1X2 annihilation; this allows the fermions X3 an d X4 after rescattering to 
still be non-relativistic and almost on-shell if their masses are close to those of the fermions in 
the initial state, thereby enhancing the loop correction. Besides, the small loop momentum q 
enables the factorization of the exchange of the boson before the annihilation, which signifi- 
cantly simplifies the calculation, as we will see later. Since in the non-relativistic approximation 
the loop correction is UV finite, no renormalization is required. 

Given the initial momenta pi,p2 and the final momenta p[,p 2 , following ref. jl] we introduce 
the four-vectors P = (pi +p 2 )/2, half the total momentum, and p = {p\ — p 2 )/2 whose spatial 
component is the three-momentum of the annihilating fermion Xi i n the cms. In this frame, P 
and p are explicitly given by, 



P = (\/iP + m 2 1 + \/ I P + ml)/2, P = 0, (3) 

and, 



Po = (y^ + mf- ^Jp 2 + m 2 2 )/2. (4) 

The momentum difference in the final state, p' = —p' 2 )/2, affects the annihilation amplitude, 
but it is irrelevant for the calculation of the correction. 



: As in ref. [5], the final result holds for bosonic WIMPs as well. 
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Figure 2: Feynman diagram for calculating the amplitude correction 5Al(\p\,p')- P = (pi + 
P2)/2, p = (pi — P2}/2, and q is the loop momentum. The coupling matrices V and f are 1 for 
scalar coupling, and the big blob represents the X3X4 annihilation vertex. 



Following the line of the argument in ref. [5], the one-loop correction term to the annihilation 
amplitude can be written as, 



SAf X2 {\p\,p) = ig^xix^lx^AP^) 



x 



(2tt) 4 [(q - P) 2 -m\ + ie][{q + P) 2 - m\ + ie] 

77 * ^ (14^) • (5) 

[{P ~ Q) ~ m i + i£\ 



Here rrii, i G {1,2,3,4} is the mass of fermion \ii an d Q is the loop momentum as illustrated in 
Figj2} The matrices T and f describe the 4>X3Xi an d 0X2X4 couplings, whose strengths are given 
by the (possibly complex) couplings g<f> XlXj - For scalar, pseudoscalar, vector, and axial vector 
couplings, T and T are (1,1), (75,75), (7^,7/J, and (7^75,7^75) (index \i is summed over), 
respectively. In the latter two cases, an extra overall minus sign should be introduced, coming 
from the propagator of the spin-1 boson 0. (75)"^ (more on this later) stands for the effective 
Lorentz structure of the annihilation process; the remaining dynamics of the annihilation is 
contained in the "reduced" amplitude A^'^dql,^). 

As in [3] the "switch" (75)^ serves to differentiate between S- [n^ = 1) and P-wave 
(til = 0) contributions, though an explanation of its validity in the extended scenario here is in 
order. In the previous work [5] the annihilating particles were assumed to be identical Majorana 
fermions. In this case the Pauli exclusion principle only allows one choice of the total spin S for 
each partial wave, and the CP parity of the initial state is fixed. Here the initial particles can 
be different fermions. The total spin S is then no longer determined uniquely by the orbital 
angular momentum L, i.e. by the partial wave considered. In general all the possibilities listed 
in Table Q] need to be considered. 

However we can opt to focus on the J = configuration as representative for a given partial 
wave, i.e. the l So state in the S"-wave and 3 Pq state in the P-wave. This can be described by 
a scalar-like effective vertex, i.e. the (75) ni Lorentz structure, including both scalar (n^ = 0) 
and pseudoscalar (n^ = 1) types, to describe the annihilation of WIMPs. Here one is exploiting 
the fact that the complicated three-momentum dependence of the Sommerfeld enhancement 
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S=0 S=l 


L = 
L = 1 


1 C 3 c 
DO Di 

1 P 1 3 P , 3 P 1 , 3 P 2 



Table 1: Possible spin states in each partial wave. For each state labeled by 2S+l Lj, S is 
the total spin, L the orbital angular momentum (partial wave), and J is the total angular 
momentum. 

depends only on the partial wave type, i.e. it is independent of S and jj^] 

A noteworthy feature of Eq.(|2]) is the mixing between different channels. 8A^ lX2 is the 
correction to the amplitude A\] X2 of the annihilation of a X1X2 pair, while A x 3X4 on the right- 
hand side is the "reduced" tree-level amplitude for X3X4 annihilation. As a result we can 
no longer hope to factorize the correction at the cross section level, unlike in ref. [5j. This 
considerably complicates the calculation of these corrections in actual applications, as we will 
demonstrate in the sample calculations of Section |4} 

We simplify Eq.g using the same approximations as in ref. [5]. First, the q dependence in 
the "reduced" bare amplitude A X3 £ 4 is neglected in the non-relativistic limit. This allows to pull 
the bare amplitude A Xi ^ 4 out of the integral, i.e. factorization still works at the amplitude level. 
The bosonic propagator, l/[(p — q) 2 — rnZ], is approximated by the instantaneous "Coulomb- 
like" part, —l/[(p — q) 2 + since in the non-relativistic limit the energy exchange is much 
smaller than the momentum exchange. The rest of the denominator has two poles beneath the 

real axis of q°, situated at U4 + P° — ie and U3 — P° — ie, where = ^jq 2 + m| 4 . The first 

pole gives a much larger denominator and its residual is neglected. The denominator of the 

residue of the second pole is: 

m 3 + m 4 n 

u = • (mi + m,2 — 1713 + m^) ■ p 

m 4 

'q 2 17131714 m\ + rri2 217131714 1 

— ■ 1 ■ -^(m 3 + i7i4-m 1 -m 2 

-P z 1713 + 1714 m\m,2 1Tl3 + 17l4p z 

x[(p-q) 2 + ml}. (6) 

Here we have performed non-relativistic expansions of all energies, keeping only the leading 
non-vanishing powers of p 2 and q 2 . Defining two auxiliary parameters, 

Co = ■ (mi + m 2 - 1713 + m 4 ), (7) 

1TI4 



x 



IIO3IIO4 IILI T 1102 *llb3H04 -L , , , . 

k = — — — — — — -(m3 + m 4 -mi-m 2 ), (8) 



17131114 rrti + ui2 271131714 1 
7n\7Jl2 7773 + 1714 7713 + 1714 p 

Eq.(|6| can be written in a succinct way, regardless of the partial wave 



V = c D -f 



p z 



q) A + mi\. (9) 



2 We will see later that the sign and strength of the Sommerfeld correction can also depend on S and J; here 
we wish to compute the loop functions describing the non-trivial dynamics of the corrections, which depend 
only on L. 
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The numerator of Eq.(j2]) differs for the two partial waves we are considering. In case of 



S-wave annihilation, 



1. In this case we can set all 3-momenta in the numerator to 



zero. Moreover, since we can get an enhanced correction only if all four participating fermions 
Xi have similar masses, we ignore terms oc {m\ — m 2 )(m 3 — m 4 ) which is of second order in 
mass differences. Performing a string of gamma matrix algebra, the Lorentz structure of the 
numerator 

N := v(p 2 )Y^ -f + m 4 )(75) ni (^ + f + m 3 )Tu( Pl 

is reduced to 

M = c^v(p 2 h 5 u(Pi) ■ 
The coefficient depends on the types of 4>XiX3 an d 



(10) 



'JV 



UX2 vertices: 

(mi + m 2 ) 2 /4 + m 3 m 4 + (m 1 + m 2 )(m 3 + m 4 )/2, 
- (mi + m 2 ) 2 /4 - m 3 m 4 + (mi + m 2 )(m 3 + m 4 )/2, 
(mi + m 2 ) 2 + 4m 3 m 4 — (m 4 + m 2 )(m 3 + m 4 ), 
-[(mi + m 2 ) 2 + 4m 3 m 4 + (mi + m 2 )(m 3 + m 4 )]. 



scalar 

pseudoscalar 
vector 
axial vector 



(12) 



Here for vector and axial vector couplings the negative sign in the propagator has been ab- 
sorbed into cfj. If we write m 2i3j4 = mi(l + 62,3,4) and expand up to linear order in the mass 
differences described by the e,, the results for the scalar and vector are the same, being given by 
4m 2 [1 + (e 2 + e 3 + e 4 ) /2]; the result for the axial vector exchange differs from that for vector 
exchange by a factor of —3 [5J, whereas for pseudoscalar interaction is of order e 2 , and can 
thus be neglected to the order we are interested in. 

The bi-spinor in TV remains finite as \p\ — > 0, as expected for an S"-wave amplitude. It can 
be combined with the reduced amplitude to give the full tree-level X3X4 annihilation amplitude. 
The correction for the S- wave X1X2 annihilation amplitude is therefore proportional to the tree 
XzXa annihilation amplitude, 



8Af^{\p\,p')\ 



l-loop 



9<t>xiX394, X2Xi c° N / m 1 m 2 
8tt 2 cd\p\ V m 3 m 4 



I s (r,K)A^(\p\,p'), 



(13) 



where the universal numerical pref actors to the numerator and denominator cfj, cd are given 
by Eq.(12) and Eq.Q, respectively. The square root in front of Is occurs because in Eq.(ll) 
the external fermions are Xi an d X2, whereas the amplitude 3 g 4 obviously refers to reactions 
with x?j an d X4 i n f ne initial state; in the relevant non-relativistic limit, the spinors simply 
reduce to the square root of the mass of the respective fermion. Finally, the function I sir, k) 
describing the dynamics of the correction is defined as 



Is(r, k) = 9fte 



x 



X — K 



In 



(1 



x] 



1 — x) 2 + r 



-dx 



(14) 



k has been defined in Eq.(p|, and r is given by 



mj 



(15) 



Note that the integral in Eq.(14) should be understood as a principal value integral. 
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Before evaluating the integral in Eq.(14), we discuss the case of P-wave annihilation, which 
corresponds to = 0. Here we again neglect terms that are of second or higher order in 
fermion mass differences, but we keep terms linear in the 3-momenta p, q. Similar algebra 
as for S'-wave annihilation shows that the numerator is proportional to v(p 2 ) (a>r4 + &r) u ipi)i 
where the constants a-p and &r depend on the Dirac structure of the 4>XiXj couplings. Note that 
both terms are of first order in the three-momentum: A is explicitly of this order, but the 7 
matrix couples the two large spinor components, as can easily be seen in Dirac representation. 
The second term is 0(p) since the bi-spinor only contains products of one large and one small 
spinor component. The term proportional to bp, which vanishes for vanishing mass splitting 
between the four fermions, can be treated straightforwardly. Since in a one-loop calculation 
we need the interference between the one-loop amplitude and the tree-level amplitude, we 
treat the term oc ar by multiplying the one-loop correction with the hermitean conjugate of 
the tree-level amplitude v(p 2 )u(pi), and dividing by the square of the tree-level amplitude. 



In other words, we replace the one-loop amplitude 5Ai t p by [SA^pA^ 
leaves the relevant product SAi^Aq p unchanged. This yields: 



A p/\A 



o,p\ 



which 



Af 



/|X3,X4 



4 



+ 4 



Q-P 

p 1 



AX3,Xi 
A 0,L > 



(16) 



where the numerical factors c^, are 



mim 2 . 

2 (mi + m 2 + m 3 + m 4 ), 

mi + m.2 

2 («3 + W4 — mi — m2j, 

mi + m 2 

„ mim 2 , , 
-(m 3 + m 4 ), 
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mi + m,2 
miW2 



k mi + 7712 



(m 3 + m 4 ) 



scalar 

pseudoscalar 
vector 
axial vector 



(17) 



a N 



' 77131114 — (mi + ^2) 2 /4, 
(mi + m2) 2 /4 — m 3 m4, 
4m 3 m 4 — (mi + m^) 2 ■, 
(mi + 7712) 2 — 4m 3 m 4 , 



scalar 

pseudoscalar 
vector 
axial vector 



The coefficients for vector and axial vector interactions again contain an extra factor of —1 from 
the sign of the spin-1 propagator. We note that the are all of the same form, but differ by 
overall factors; they all vanish linearly for vanishing mass differences. To linear order in mass 
differences, for vector and axial vector interactions is the same as cf r for scalar or vector 



N 



interactions, but for scalar interactions differs, and cj^ for pseudoscalar interactions vanishes 
only linearly in mass differences. A very light pseudoscalar with off-diagonal couplings^] could 
therefore give significant corrections to co-annihilation. 



For pseudoscalar coupling, c 
777,4 = fni- 



^ and 



vanish if the coupling is diagonal, in which case 777.3 



777! and 
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In the end, the P-wave amplitude correction can be written as 



(19) 

The square root factor occurs for the same reason as in eq.(13), and the function Ip(r, k) is 
defined as 



Ip(r, k) = 5fe 



2x z 



X — K 



1 + 



x 2 + 1 + r , (x + l) 2 + r 
m 



Ax 



x — l) 2 + r. 



ci.r 



(20) 



So far we have assumed that X3X4 annihilation proceeds through a (pseudo) scalar vertex. 
This describes annihilation from a state with total angular momentum J = 0. As noted 
earlier, if the initial state consists of two identical Majorana fermions, there is a one-to-one 
correspondence between orbital angular momentum L and spin S, such that L = (S'-wave) 
requires spin S = 0, and hence J = 0. However, we saw above that this need no longer be 
true for co-annihilation processes. We therefore repeated the S'-wave calculation for the case 
that X3X4 annihilate through a 7^ vertex. This describes annihilation from a J = 1 state. 
Only a space-like index, v = k G {1,2,3}, gives a non-vanishing tree-level amplitude in the 
limit of vanishing three-momentum. We find that the coefficients cfj for scalar, vector or axial 
vector interaction of the exchanged boson are now all equal to for scalar boson exchange and 



annihilation through a 75 vertex, as given in Eq.(12); the coefficient for pseudoscalar interaction 



again vanishes, up to terms that are quadratic in mass splittings. In particular, in this case no 
factor —3 appears for axial vector exchange. Note that L = and J = 1 implies S = 1. This 
is consistent with the rescattering argument of ref . [5 J . 





scalar 


pseudoscalar 


vector 


axial vector 


spin singlet 


1 





1 


-3 


spin triplet 


1 





1 


1 



Table 2: cs or ctv/(cd^ x ) ^ or different spin states of the initial fermion pairs with different 
boson-WIMP coupling types. 

In the limit of vanishing mass splitting we can therefore write the overall factor cn / '{cdtti^ 
for either the S- or P-wave as shown in Table [2j This coefficient simply becomes 1 (0 for 
pseudoscalar exchange), except for the case of axial vector boson exchange in a spin-singlet 
(S = 0) state. We see that this coefficient is in fact determined by S, rather than by L and 
J. Again in the limit where the fermion mass splittings can be neglected everywhere except 
in the definition of k appearing in the loop functions, the correction to the X1X2 annihilation 
amplitude can then be written as: 

SAf * 2 (b1,p')li-Loo P = cs 9 *™^ ^hir, k)A$?(\p\,p') , (21) 

where m x is the mass of the WIMP. Recall that we are interested in the calculation of pertur- 
bative one-loop corrections. We expect the overall magnitude of these corrections to be of the 
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order of at most 10 or 20%; for larger corrections, resummations will be necessary. Moreover, 
co-annihilation is important only for mass splittings below 10% or so; in fact, in our numerical 
examples we will encounter much smaller mass splittings, as relevant for the annihilation of 
higgsino- or wino-like states in the MSSM. Contributions of order {a/ii){5m/m x ) can then 
safely be neglected. On the other hand, Eq.(|8]) shows that the quantity k appearing in the 
denominators of the loop functions will diverge for any finite mass splitting when the initial 
three-momentum p — > 0. It is therefore important to take the mass splitting into account when 
computing k. 

Note that Eq.(21) is applicable also to the (co-) annihilation of bosonic WIMPs, with spin 
or 1. This has been shown in [5] for WIMP self-annihilation, and remains true also for the 
more complicated situation discussed here. 

We have found analytical expressions for the integrals iz(r, k) using contour integral meth- 
ods, as follows. In the S partial wave, 



Is(r, k) 



C s {r, k), 

C s (r, K) + P s (r,K) 



K>0, 
K < 0, 



where C*s(r, k) comes from the branch cut of the logarithm in Eq.(14): 

2^/r 



% ■ arctan 



C s (r,K) 



( 2 Vr \ 
ix ■ (arctan ^ — - — J + ir), k < — r + 1, 



,k — 1 + r) ' 
2^ 
+ 

and Ps(r, k) is the residual at the pole iyJ—K when k is negative 

2 



k > — r + 1, 



P s (r, k) 



In the P partial wave, 



7r • arctan 
7r ■ (arctan 



-K 



K + 1 + r 
2 v /3 k 
k + 1 + r 



7T 



K > —r — 1, 
K < — r — 1. 



C P (r, k), 

C P (r, K) + Pp(r,K) 



K > 0, 
K < 0, 



(22) 



(23) 



(24) 



(25) 



where Cp(r, k) comes from the branch cut of the logarithm in Eq.(20): 

2^ 



C P (r, k) 



71 



r + 



k + 1 + r 



arctan 

K — 1 + r 

arctan : h tt, 



k — 1 + r 

and Pp(r, k) is the residual at the pole iyf—K when k is negative 

2J^k 



P P (r, k) 



-7T 



k + 1 + r 



arctan 



-K 



k + 1 + r ' 



arctan 



-K 



k + 1 + r 



+ 7T. 



K > —r + 1 ^ 
K < — r + 1 



(26) 



k > —r — \ 
K < — r — 1 



(27) 
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We note that the "classical" Sommerfeld effect refers to the exchange of a massless boson 
(i.e., r = 0) between fermions of equal mass (i.e., k = 1). In this case one simply has ls(0, 1) = 
lp(0, 1) = 7r 2 /2. However, for r ^ the corrections to S- and P-wave annihilation differ 
significantly, as we will see shortly. 

2.1 Complications due to Fermion Flow and Spin 

Our discussion so far assumed implicitly that "Dirac arrows" can be drawn consistently along 
the fermion line, allowing to directly read off the correct order of external spinors, propagators 
and vertex factors. This is always true in the SM, thanks to the "accidental" conservation 
of lepton and baryon numbers, but need not be true in extensions of the SM. In particular, 
"clashing arrows" frequently occur in supersymmetric extensions of the SM [10J. 

In our case these occur in particular in diagrams with Majorana fermions; an example is 
shown in Fig. [3j We use the convention of Denner [TTJ to systematically treat such Feynman 
diagrams. In this treatment one introduces an auxiliary fermion flow, which is continuous 
through the diagram, as shown in the right diagram of Fig. [3j This auxiliary fermion flow is 
used to write down the spinor chain for this diagram. In most cases the original vertex factors 
should be used; however, if the auxiliary fermion flow goes against the direction of the usual 
Dirac arrow on a given vector-fermion-fermion vertex, with Dirac structure 7^, then this vertex 
receives an additional minus signj^] If this procedure changes the order of external spinors, the 
whole amplitude should be multiplied with — 1. This rule ensures that the direction one chooses 
for the auxiliary fermion flow is not relevant. Moreover, since the only modification required is 
a possible minus sign in front of the amplitude, the calculations in the main part of this Section 
are not affected. More details on this method, and several examples, can be found in ref.|llj. 




Figure 3: At the left is an example of a relevant Feynman diagram with clashing Dirac arrows, 
describing a one-loop correction to the s-channel annihilation of a chargino and a neutralino 
in a supersymmetric extension of the SM; here xt nas been defined as "particle" , with Dirac 
arrow along the momentum direction. At the right is the same diagram with an (arbitrarily 
chosen) auxiliary, continuous fermion flow. 

A second complication occurs when the products 9<t>xiX39<i>X2Xi anc ^ 9<t>x\x^<l>xiX3 are b°th 
nonzero and X3 7^ Xi- I n such cases there are two contributing Feynman diagrams where the 
intermediate particles are swapped; a pair of examples is shown in Fig. |4| In this case one needs 
to distinguish between X3X4 an d X4X3 annihilation. Of course, at the level of cross sections these 

4 No such extra sign appears for axial vector vertices, with Dirac structure 7^75; the difference is due to the 
different behavior of these two Dirac structures when sandwiched with the appropriate combination of charge 
conjugation matrices [llj . 



10 



w + 



YVwwvv/ 




-T^ WWVVV n- 



w + 



— — fwwwv 




Figure 4: The annihilation of X1X1 with the intermediate state XiXi can have two contributing 
Feynman diagrams. Again xt is a "particle", i.e. Xi is an "antiparticle" and has Dirac arrow 
opposite to the momentum direction. The auxiliary fermion flows are added in accordance with 
the Denner's convention. 



are the same (if consistent definitions of the scattering angle are used in both cases), but the 
corresponding amplitudes may differ by a sign. This sign can be determined as follows. 

Let A(p, cos 9, . . . ) be the reduced amplitude for X3Xi annihilation; here p is the absolute 
value of the cms three-momentum in the initial state, 9 is the cms scattering angle, and 
. . . denotes possible other quantum numbers (e.g. the spin). The reduced amplitude of the 
"crossed" diagram, for XiXz annihilation, is then given by SA(p, — cos 9, . . . ). The sign of cos 9 
has to be changed since by convention the first annihilating particle has a fixed direction: if 
X3 goes in +z direction in X3Xi annihilation, it goes in —z direction in XiX3 annihilation. The 
overall sign S depends on the spins involved. Note first of all that a crossed intermediate state 
only appears if all four Xi obey the same statistics (Bose-Einstein or Fermi-Dirac) Q If they 
are fermionic, S contains one factor of —1 from the crossing of fermion lines. An additional 
factor arises from the symmetry of the spin wave function]^ if Xz an d Xi are both spin— 1/2 or 
spin— 1 particles. In the former case the spin wave function is symmetric for total spin 5 = 1 
and antisymmetric for 5 = 0. If both X3 an d Xi have spin 1, then the spin wave function is 
symmetric for S = or 2, and antisymmetric for S — 1. Either way, an antisymmetric spin wave 
function leads to an additional —1 factor in S. As a result, S = (— l) s for the annihilation of 
either two spin— 1/2 fermions or two spin— 1 bosons; however, the co-annihilation of one scalar 
and one vector boson always gives S = +1. 

2.2 Fermion Exchange 

Before concluding this Section, we briefly comment on the exchange of light fermions between 
co-annihilating WIMPs. This is possible, e.g., in the co-annihilation of a supersymmetric 
neutralino and a slepton [12J, where the exchanged fermion could be very light (a charged 

5 We will show in the next Subsection that the particle exchanged in the X1X2 — * X3Xi rescattering has to 
be a boson. 

6 The symmetry of the orbital angular momentum part of the wave function has already been described by 
the change of sign of cos 9. 
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fermion) or even nearly massless (a neutrino). Moreover, since the "Sommerfeld correction" is 
non-relativistic in nature, one might naively expect that exchange of a light fermion also leads 
to enhanced corrections; after all, the spin of the annihilating WIMPs does not matter in the 
usual Sommerfeld corrections. 

However, one can see fairly easily that the exchange of a light fermion between (co- 
annihilating WIMPs does not yield an enhanced correction. The reason is that the numerator 
of the fermion propagator contributes an extra factor — $ + rrif to the numerator of the ar- 
gument of the loop integral. Recall that q has to be counted as being of order \p\ here, and 
qo oc \p\ 2 . If the fermion is light, all three terms in this extra factor are therefore or 
smaller; this implies that the correction will not be enhanced by a 1/v factor even for rrif — > 0. 
On the other hand, if mj is comparable to the WIMP mass, this extra factor in the numerator 
does not give a significant suppression; however, in that case the loop integral is small, as 
general arguments in the Introduction indicate and the numerical results in the next Section 
show explicitly. Hence there is no range of mj where one can expect an enhanced correction 
from fermion exchange. We therefore do not discuss these contributions any further. 



3 Discussion of Results 



The solutions Eq.(22) and Eq.(25) are two model-independent functions of variables r and k 
only. In order to better understand the dependence of the size of the correction on r and k, we 
define the amplitude enhancement function 

£ L (r,n) := -L^I L ( r ,K), (28) 
This allows to recast the the amplitude correction formula as: 

<^ 2 (b1,p')|i-ioop = Cs ^£ L (r,K)A^T(\p\,p') } (29) 

where the relative enhancement of the annihilation amplitude simply consists of three parts: the 
uniform prefactor cs given in Table[2j the one-loop factor am x /m ( j } with a = 94>x\xz3%x2xJ (^ 7r )' 
and the amplitude enhancement function £l(v, k). Here and in the following numerical results 
we work to zeroth order in WIMP mass differences wherever possible, i.e. we set mi = m<i = 
ms = m4 = m x everywhere except in the definition of k, Eq.(|8]). 

We first discuss several examples in order to illustrate the size of the radiative corrections 
we are calculating, and to understand the physics. We take 80 GeV (the mass of the H^-boson) 
for the mass of the exchanged boson, and 1.1 TeV for the mass of the lightest among the four 
fermions (i.e. the dark matter particle); the latter is roughly the mass required by the thermal 
relic density today if the WIMP is a higgsino-like neutralino. Moreover, for simplicity we 
consider only two co-annihilating states. Since we are interested in situations with small mass 
splitting, m 2 — vn\ <C m 1; we essentially have to consider only three processes. 

The first reaction, xi + Xi ~ > Xi + Xi — X + Y, stands for all reactions where the fermion 
masses before and after rescattering are the same, i.e. these results are applicable (with very 
small changes) to any process of the kind XiXj ~^ XiXj — ^ X + Y. This is the case discussed 
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in [5]. Here k = lQso that 



£s(r, K — 1) = ^— ■ arctan (—t= \ , (30) 
2 Vyr/ 

£ P (r, k = 1) = ^- [- v^+(1 + r/2) arctan (-^1 . (31) 
2 L V,/r/J 



The Sl(t, k = 1) are plotted as functions of \p\ for both the 5- and P-wave in Fig. 5a We see 
that £s saturates at 1 and Sp saturates at 1/3 in the zero-velocity limit, \p\ — > 0. For larger 
\p\ Ss remains larger than Sp, but the two functions approach each other for large \p\. Recall 
that \p\ ^> corresponds to r <C 1, where Is(r, k = 1) and Ip{r, k = 1) both approach 7r 2 /2, 
i.e. in this limit Ss(r <C 1, K = 1) = £p(r <1,k=1)4 7rm0/(4|p|). We finally note that Sp 
has a very broad maximum at \p\ ~ m^/2; however, the value at this maximum exceeds the 
value for \p\ — > only by 6.8%. Nevertheless this maximum at non-vanishing \p\ implies that 
Sp remains approximately constant out to much larger momenta \p\ than Ss doesj^] 

Next we consider the process where the particles in the intermediate states are heavier than 

the initial ones: Xl + Xl ~ * X2 + X2 ■> X + Y, with m 2 = mi + 5m. Since our correction 

function depend primarily on the total mass difference rrii + rrij — — mi , one finds very similar 
results for X1X1 ~^ X1X2 if the mass difference 5m is doubled. 

This case differs from the one we just discussed since now 1 and is no longer a constant. 
Note that the entire dependence of the correction functions on the mass splitting is described 
by this parameter. For small mass splitting, \5m\ m\, k can be written as 

5m ( 2m? \ , , 

« - 1 + — 1 - > 32 

m\ \ V J 

where 25m is the difference between the sum of the masses in the initial state and the sum of 
the masses in the intermediate state; in the case at hand, 5m = m 2 — m±. In the limit p — > 0, 
k can be further simplified to 

5m 2m? 

« ^ ^t- 33 

mi p z 

Note that for the reaction we are discussing, k — > — 00 as \p\ — > 0. Therefore even very small 
mass splittings have to be kept, if we want to describe the radiative corrections correctly at all 
values of 



The corresponding amplitude enhancement factors are plotted in Fig. 5b Comparison with 
the first case discussed above shows that the correction still reaches a plateau as \p\ — > 0, albeit 
at a reduced value. This can be understood by expanding the functions Sl{t,k) in terms of 
5m 1 'm in the limit \p\ — > 0. In the scenario we are considering, both the functions Cl given in 



Eqs.(23) and (26) and the functions Pl given in Eqs.(24) and (27) have to be included. This 
gives: 

£s(r,/0l| Pl ^ - !=^r> ( 34 ) 



7 In scenarios with three of more co-annihilating WIMPs it is possible to have rescatterings XiXj ~> XkXi 
such that rrii + m,j = + mi but rriimj ^ m^m;, in which case k is still independent of \p\ but differs from 1. 
8 This maximum is not reproduced by the numerical approximation of ref . [S] . 
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(c) X2 + X2 -> Xi + Xi -> X + Y 

Figure 5: The amplitude enhancement factors Si(r, k) are plotted as a function of \p\ in the CM 
frame for the cases that the intermediate state is the same as the initial state (top), is heavier 
than the initial state (middle), and is lighter than the initial state (bottom); in the latter two 
cases, results for two different mass splittings are shown. These results are for a WIMP mass 
of 1.1 TeV and a boson mass of 80 GeV. 
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Note that 2m\5m > m\ for all cases considered, leading to a sizable suppression of the correc- 
tion, especially for S'-wave annihilation. The suppression is less in the P-wave case because of 
the extra factor in parentheses. 



Another characteristic feature of the curves in Fig. 5b is the occurrence of pronounced 
maxima at the threshold value of \p\ where the intermediate state can be produced on-shell 
in non-relativistic kinematics. This happens at the point \p\ 2 = 2mi5m, which corresponds to 
K — 0. The S-wave function £$ has a cusp at this point, i.e. is not different iable, while Sp as 
a function remains smooth at this maximum. 

The physics is therefore clear. When \p\ is below the threshold for real \2 pair production, 
the intermediate particles \2 are produced virtually and the correction is suppressed. At the 
threshold the intermediate state can finally be produced on-shell with zero relative velocity 
and the propagators of \i m the one-loop correction are large. Afterwards the correction 
decreases again with increasing momentum. Note that for \p\ values near the maximum, the 
loop correction for a heavier intermediate state can exceed that for diagonal scattering, i.e. for 
k = 1, discussed above. 

We finally discuss the case where the particles in the intermediate state are lighter than 
those in the initial state, i.e. X2 + Xi — > Xi + Xi — X + Y . Again, results for X2X2 — > X1X2 
are very similar, if the mass difference 5m is increased by a factor of 2. The resulting loop 



functions £l(\p\) are plotted in Fig. 5c 



We again observe plateaus for \p\ — > 0. The finite mass splitting again leads to a suppression 
of the correction functions in this limit. Note that now k — > +00 in this limit, so that the Pl 



functions of Eqs.(24) and (27) do not contribute. This leads to a stronger suppression than 



in the previous case where the particles in the loop were heavier than the external particles. 
Expanding the correction functions in the mass splitting, which is now negative, we find: 



£s{r, k)\ 



1 



, mi \Sm\ 



£ P (r, k)\ 



'1 + 3-2 



mi \Sm\ 

2 



7)1 



4> 



Ip1->o ~ 3 ]i 2 mi \ 5m \ \ 1 4- 2 mi l <?m l 

4> 



(36) 



(37) 



The enhancement function £$ is suppressed more strongly than for positive dm, once 2m\\5m\ > 
m 2 ,. The suppression for the P-wave is again weaker than for the S- wave. 



Note that the expansions (34), (35), (36) and (37) reproduce the exact corrections rather 



well as long as the kinetic energy (not the momentum) in the initial state is smaller than the 
absolute value of the mass splitting, i.e. for p 2 < 2mi\5m\. For the larger mass splitting shown 
in Figs. 5b and 5c 5m = 100 GeV, this remains true for nearly the entire momentum range 
where the non-relativistic expansion can be trusted. 

In summary, when the mass splitting is not vanishing, the correction for small initial mo- 
menta is suppressed. This suppression is stronger if the particles in the loop are lighter than 
the external particles, and always increases with the absolute value of the mass splitting. To 
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Figure 6: Contour plots of the amplitude enhancement function El in the (k, r) plane for S'-wave 
(left) and P-wave (right). 



illustrate this, we calculate the value of the mass splitting \5m\ where the correction for \p\ — > 
is suppressed to 10% of the correction for 5m = 0. Using formulae (34) and (36) for the S'-wave, 
the corresponding relative mass splitting \5m\/mi is about 20% and 2% for positive 5m and 
negative 5m respectively. Intermediate states with yet larger mass splitting can be deemed as 
relatively unimportant for the correction to the annihilation cross section. 

Before concluding this Section, we give contour plots of the amplitude enhancement function 
El in the (k, r)— plane, for both S and P partial waves (Fig. [6]). In both plots, El is large where 
r is large and the absolute value of k is small, corresponding to small three-momentum and 
small mass splitting. The maximum value of £l is 1 and 1/3 for S- and P-wave, respectively. 
Recall that the magnitude of the enhancement also depends on the one-loop factor am x /m ( f ) . 
In order to produce a sizeable correction, m x /m ( f ) should be large enough while at the same 
time not too large, in order to keep am x jm$ below 1 so that the one-loop approximation makes 
sense. In the framework of WIMPs, where the coupling constant is weak, there is still a fairly 
large range of mass for dark matter particle satisfying this condition. 



4 MSSM in the Co— annihilation Region 

The results presented so far are model independent. However, in order to gauge the importance 
of the corrections to the annihilation rates and the current relic density, we need to perform 
numerical computations in the framework of a specific WIMP model where some particles have 
masses close to that of the dark matter particle and therefore co-annihilation can happen. In 
this chapter we choose the Minimal Supersymmetric extension of the Standard Model (MSSM) 
for this purpose. Before discussing one-loop corrections, we give a brief review of the properties 
of neutralinos and charginos, as well as the form of the Boltzmann equation we will work with. 
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We then describe the numerical implementation of the corrections, before presenting numerical 
results. 



4.1 Formalism 

In the MSSM the four neutralinos are linear combinations of four different neutral fermionic 
interaction eigenstates: the bino B, the neutral wino W°, and two neutral higgsinos H±, H®. 
These states mix as a result of electroweak gauge symmetry breaking. This mixing can be 
described by the mass matrix of the four fermions in the basis (B, W°, Hf, H®) [10J, 



/Mi -M z s w c(3 M z s w sf] \ 

M 2 M z c w c p -M z c w s p 

-M z s w c p M z c w C(3 -fi 

\ M z s w sp -M z c w s p -ft J 



(3f 



Here, M\ and M 2 are the bino and wino masses, respectively, and \i is the supersymmetric 
higgsino mass parameter. The off-diagonal terms, which cause higgsino-gaugino mixing, result 
from gauge-strength Higgs-higgsino-gaugino interactions, which contribute to the neutralino 
mass matrix when the Higgs fields attain vacuum expectation values (vevs). sw is shorthand 
for sin 9w, and cw, sp, cp stand for cos 6^, sin/3, and cos/3. Here tan/3 is the ratio of the vevs 
of the two Higgs fields. The mass matrix Ai° is diagonalized by a 4 x 4 unitary matrix Z to 
produce four Majorana neutralino mass eigenstates, 

X° = Z a B + Z i2 W° + Z a Hl + Z l4 H° 2 . i = 1, 2, 3, 4 (39) 

Which one of the four neutralinos is the lightest depends on the parameters M 1; M 2 , /x and tan /3. 
The lightest neutralino can be a good WIMP Dark Matter candidate, if it is stable. This 
condition is satisfied if the lightest neutralino is also the lightest of all superparticles (LSP) 
and if R— parity, or a similar symmetry under which particles and superpartners transform 
differently, is preserved. 



The parameters that determine the neutralino mass matrix (38) also appear in the Dirac 



mass matrix mixing the charged wino and higgsino states. It is given by [TU] 

M ± = ( r M * ^ M ^\ (40) 

Since the chargino mass matrix is not symmetric, one needs two unitary matrices U, V for its 
diagonalization, i.e. the left- and right-handed components of the chargino mass eigenstates 
X12 m i x differently [TU]. 

Co-annihilation between the lightest neutralino and lightest chargino is important when- 
eveiQ M 2 < |Mi|, \fi\ or < \Mi\, M 2 . In the former case both and xt are dominated 
by their wino components. To good approximation these three states form a triplet under the 
weak SU(2) gauge symmetry. The tree-level mass difference between the lightest neutralino 
and chargino is then of order M z /Mifi 2 [T3]. This is extremely small (Mi, fi ^> M z ), and radia- 
tive corrections have to be included. If scalars are somewhat heavier than the lighter wino-like 



9 Wc follow the usual convention where Mi is real and positive; this can be assured by phase transformations 
without loss of generality. M\ and/or /i can then be negative or, in the presence of CP-violation, complex. 
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states, these corrections amount to about 170 MeV [IB] , so that the chargino-neutralino mass 
differences is almost independent of the parameters of the neutralino mass matrix. 

In the second scenario, Xi> Xi an d X2 are au dominated by their higgsino components. Up 
to small corrections, they form two two-component doublets of the SU(2) gauge symmetry, 
which can be grouped into a single doublet of Dirac fermions. The mass differences are of order 
of M|/Mi, MI/M2 [H], which are still small, but much larger than in the wino-dominated 
case. In the presence of large ti — £r mixing radiative corrections to the mass splitting can be 
of comparable size as the tree-level splitting |I5[[T5]. However, in the region of parameter space 
where the thermal relic density of a higgsino-like Xi has the right magnitude, the relative mass 
splittings remain quite small even after one-loop corrections. Here we ignore these corrections 
for simplicity. Note also that in this scenario the mass of the lighter chargino is typically about 
midway between the masses of \i an d X2- 

In both cases, other fermionic particles exist whose masses are close to that of the dark 
matter particle. The coannihilation mechanism becomes important, including the off-diagonal 
Sommerfeld effect studied in this paper. Recall that the latter effect is very sensitive to the 
mass splitting, as observed in the previous Section. For wino-like WIMP the mass splitting is 
so small that it can be neglected in the calculation of the relic density. On the other hand, for 
higgsino-like states the mass splitting can be significant, although the relative mass splitting 
Sm/\n\ decreases rather quickly with increasing WIMP mass ~ 

Another special phenomenon in the case of higgsino dominance is the physical phase be- 
tween Xi an d X.2- Simply diagonalizingp^j the bottom right 2x2 block (the higgsino sector) of 
the neutralino matrix one finds that one of the two mass eigenvalues becomes negative, if 
Zij is kept real. One can then multiply the fermion field with apparently negative mass with a 
factor «7 5 , leading to a state with positive (i.e., physical) mass. Then the Feynman rules that 
involve this particular fermion field need to be modified accordingly [To] . It can be cumbersome 
to keep track of this special fermion field in calculations. The other way is to relinquish the 
reality constraint on the matrix and multiply the row of associated with the field with 
otherwise negative mass by an imaginary unit %. We choose to adopt this convention. As the 
elements of the mixing matrix Z^ appear ubiquitously in Feynman rules for vertices involving 
neutralinos, the physical relative phase can have a significant impact on the Sommerfeld calcu- 
lation, as we will see later. This phenomenon does not affect the chargino sector, since the two 
diagonalization matrices U and V allow enough freedom to make all chargino masses positive 
even if U and V are real, as long as CP is conserved in the chargino sector. 

The rest of this Subsection gives a brief review of the formalism of coannihilation calculation 
that we use, closely following refs.jH] and [T7j, where further details can be found. 

Consider a chain of supersymmetric particles Xi = 1, iV) whose masses are close: 
mi ^ m-2 ^ ■ ■ ■ ^ rriN-i ^ n^N (xi is the WIMP). Then the scatterings 

Xi + X^xj + Y (41) 

will be frequent enough to maintain the relative equilibrium between the densities of these 
particles even long after they have collectively decoupled from the thermal bath of the standard 
model particles, provided the mass splitting between the heaviest and the lightest of these 

10 The diagonalization of the matrix M a is Mq = Z*M°Z~ 1 ; see for example Section 9.2 in [TU] , 
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sparticles is significantly smaller than m\. This implies 



eq 

^ = ^ . (42) 

n n C( i v ' 

However, at some time after the decoupling of xi & h heavier states will decay into Xit which 
is our dark matter candidate particle. We therefore only need to keep track of the sum of 
all superparticle densities, n = ^n^. Consequently the Boltzmann equation describing the 
evolution of the number density of the dark matter particle is augmented from, 

= -ZHn^ - (a m v) (n| - n e f) (43) 

to, 

^ = -3Hn - K%) ( n Xi n Xj - ) ' ( 44 ) 

where oij is the cross section for the annihilation of Xi an d Xj m t° Standard Model particles and 
(...) denotes thermal averaging. Here we have assumed that all Xi remain in kinetic equilibrium 
during the epoch of chemical decoupling; this is usually the case, since elastic scattering of Xi 
particles on SM particles are much more frequent than XiXj annihilation reactions. 

The effects of all the coannihilation channels XiXj ~ > XY can be encapsulated in a new 
quantity, the effective cross section a e s: 

n eq n eq A 

M = J2(^)^^-. (45) 



ij 



eq 



This allows to recast the Boltzmann equation in a succinct way similar to the expression (43) 
without coannihilation: 

^ = -3Hn - (a cS v) (n 2 - n 2 cq ) . (46) 

Basically, (a e sv ) is just a weighted sum of cross sections of many (co-) annihilation processes. 
The equilibrium total number density n eq in the denominator of the last expression in Eq.(45) 
is, using the Maxwell-Boltzmann distribution for ff. 



n eq 



T J2 9i mlK 2 {^), (47) 



2tt 2 



where gi is the number of internal degrees of freedom of xi- Similarly, the numerator A can 
simplified to, 



A =%f dp^W^K, ) . (48) 



The functions Ki(x), K2(x) appearing in Eqs.(47) and (48) are the modified Bessel function of 



the second kind of order one and two respectively, and p e s is the absolute value of the three- 
momentum of xi m the center-of-mass frame of the XiXi pair, so that s = 4(m 2 +p 2 s )- Finally, 
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W e g is the dimensionless effective annihilation rate that contains weighted contributions from 
every (co-) annihilation channel: 



W, 



eff 




{raj - rrij) 2 } [s - (m 8 + mj ) 2 } g igj 

,2 *i - 



s(s — 4mf) 
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(49) 



The dimensionless (unpolarized) annihilation rate per unit volume Wij is normalized to 2Ei-2Ej 
and is related to the (unpolarized) cross section via 



Wi. 



AEiEjCTijVij. 



(50) 



The square root factor in eq.(49) is understood to imply that the contribution from vanishes 
if \fs < mi + rrij . 

The remaining task is thus the calculation of the a^. These define W e g via Eq.(49), which 
in turn allows to calculate (cr e gv) via Eqs.(45)-(48). The modified Boltzmann equation (46) 
can then be integrated numerically. This procedure is the basis of numerical packages like 
micrOMEGAs dHHH] and DarkSUSY HUE]. 



4.2 Numerical Implementation 

We use DarkSUSY [20, 21J to implement our corrections to the (co-) annihilation cross sections, 
and hence to the predicted x1 Te ^ c density, since it provides separate subroutines (in FORTRAN) 
for the calculation of all helicity amplitudes for any annihilation process one cares to include 
in the analysis. We saw in Sec. 2 that the Sommerfeld corrections to co-annihilation only 
factorize on the amplitude level, not on the cross section level, so we need all relevant amplitudes 
including their phases. Moreover, the discussion of Table 2 showed that we have to keep track 
of the spins in the initial state; similarly, we saw in Sec. 2.1 that in some cases the sign of 
interfering amplitudes depends on the spin, or total angular momentum, of the intermediate 
state. 



We see from Eqs.(13) and (19), keeping the finite mass splitting between co-annihilating 
neutralinos and charginos only in the coefficient k defined in Eq.(|8]), that the one-loop corrected 
(co-) annihilation amplitude into a given final state consisting of two SM particles can be written 

as 



am^o 



Ai(b1,p')li-ioop = Ai iL (\p\,p') + £ CjV — {K(i,j),r(mt)) 4 L (b1,p') . (51) 

Here the index i (j) labels the initial (intermediate) state consisting of two x fermions, and A l 
(Aq) is the corresponding tree-level (co-) annihilation amplitude. For a particular initial state 
i, often more than one kind of intermediate state j can contribute to the one-loop correction 
via the exchange of some boson (p. Moreover, the same intermediate state j might be accessible 
through the exchange of several different (relatively) light bosons <fi. The contributions from all 



possible intermediate states and all possible exchanged bosons should be summed up in Eq.(51 ) 
to account for the complete one-loop correction. 

Let us illustrate this with a couple of concrete examples. First, consider the annihilation 
reaction x^ + Xj — ^ W + + W~ in the (more complicated) case where the (co-) annihilating states 
are higgsino-like. This is in fact one of the most important final states. In this scenario we 
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have to consider all combinations of % and j with i,j G {1, 2}. The initial neutralino pair x^Xj 
can "rescatter" into XmXn y i a ^ ne exchange of a Z or neutral (CP-even) Higgs boson, where in 
principle again all combinations m, n G {1,2} have to be taken into account^ or it can change 
into XiXi Yl& the exchange of a W or charged Higgs boson. In either case the intermediate 
XmXn or XiXi state then annihilates into a W + W~ pair. 

As a second example, consider XiXi ~^ ZW ± (i G {1, 2}) , which is one of the dominant co- 
annihilation reactions. Charge conservation implies that only XmXi ( m {1)2}) intermediate 
states can contribute, but these states are accessible both through the exchange of a neutral 
gauge or Higgs boson coupling Xi to and through the exchange of a charged gauge or Higgs 
boson coupling x® to xf- 



We square Eq.(51) and sum up the helicities to get the one-loop corrected squared ampli- 



tude, which is proportional to the differential cross section: 

£ K(b1,p')lLoo P = £ KJ 2 + («& jV(»*)) 5> e (^r 1 ^^) • (52) 

hh j,4> hh ^ ^ ' 

Here /i and /i are the helicities of the initial particles. This expression shows explicitly that we 
need the full amplitude information, including all (relative) phases between different amplitudes, 



in order to calculate the corrections. We insert the result of Eq. ( 52 ) back into the subroutine 
in DarkSUSY that computes the relic density. 

Note that DarkSUSY does not expand the (co-) annihilation cross sections in powers of the 
initial three-momentum or, equivalently, into partial waves. On the one hand, this allows us 
to immediately use the (numerical) subroutines of DarkSUSY for the calculation of the relic 
density from the one-loop corrected annihilation cross section; recall that the corrected cross 
section cannot be cast into the usual form a = a + bv 2 . 

On the other hand, we saw in Sec. 3 that, as in the case without co-annihilation [31 21 
the Sommerfeld corrections differ significantly for S- and P-wave annihilation. For the 
purpose of computing the correction, we therefore do decompose the amplitudes into S- and 
P-wave terms by invoking the subroutine that calculates a given helicity amplitude twice, 
the first time with zero momentum (i.e. for annihilation at rest), the second time with the 
actual three-momentum in question. The first call obviously gives the constant (momentum- 
independent) contribution to this amplitude, which we equate with the S'-wave contribution; to 
the accuracy of our calculation in Sec. 2, where we only kept the leading (necessary) powers of 
initial three-momentum, this identification is exact. The entire three-momentum dependence of 
the amplitude is then assumed to be from the P-wave contribution, i.e. we assume the amplitude 
to be linear in the three-momentum when extracting the P-wave contribution by subtracting 
the result of the first call from that of the second call of the subroutine. This is not quite correct. 
In general the amplitude will also contain S'-wave contributions that depend quadratically on 
the three-momentum. Our extraction of the P-wave contribution to a given amplitude will 
therefore be correct only if the S'-wave term is suppressed, or, for roughly comparable S- and 



11 In practice the exchange of Higgs bosons yields very small corrections in the scenarios we consider, since 
Higgs bosons couple to \ states only via higgsino-gaugino mixing. Moreover, among the ZxiX^j couplings 
in the case at hand only the off-diagonal Z\iX2 coupling is sizable. In practice there is therefore only one 
combination of m, n that contributes significantly for each given combination However, our numerical 
analysis also includes all sub-leading contributions. 
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P-wave contributions, for sufficiently small three-momenta. Fortunately these are precisely the 
two cases where the P-wave contribution can be expected to be significant. In order to improve 
on this approximation, one would also have to allow additional factors of three-momentum in 
the calculation of the loop functions £l, which would add further complications without great 
improvement of accuracy. Note finally that we need this decomposition into S- and P-wave 
only for deciding which of the two loop functions is applicable; otherwise the exact momentum 
dependence of the amplitudes provided by DarkSUSY is kept. 

At this point a warning to users of DarkSUSY might be in order. While performing the 
numerical calculations described in the following Subsection, we noticed that the predictions 
of DarkSUSY for the annihilation rates of several channels, including important reactions like 
Xi + Xi —> Z + Z and Xi + Xi ~ * W + + W~ , violated unitarity quite badly for large LSP mass. 
This is illustrated by the black (solid) curve in Fig. [7J which shows va(xiXi ~~ >" W + + W~) as 
a function of the LSP mass for very small initial three-momentum \p\ = 10~ 3 m^o. At such a 
small value of \p\ basically only S'-wave annihilation contributes. Unitarity dictates that well 
above all thresholds, the cross section for a fixed partial rate should decrease like 1/s, i.e. like 
l/m| . Instead the original DarkSUSY predicted a cross section that fell for m^o < 1 TeV, but 

then started to rise again J 12 ] 

We located the source of the problem to be the unnecessary and erroneous introduction 
of imaginary parts to the t— and w-channel propagators in the expressions for the relevant 
amplitudes. Note that in the final states containing two massive gauge bosons the momentum 
exchanged in the t— or u-channel is always space-like, hence these propagators do not have 
absorptive parts. After removing these imaginary parts the cross section shows the expected 
scaling with LSP mass, as shown by the dashed red line in Fig. fff^\ Of course, we use this 
modified version of DarkSUSY for the calculation of the one-loop corrections. 

Finally, we note that in order to make sure that the one-loop correction is perturbative, 
we only consider scenarios where 2am^o/m^ < 1. We estimate the resulting upper bound on 
m^o by using the weak coupling constant aw = a era (M z )/ sin 2 6 W ~ 0.034 and assuming the 
mediating boson to be W , which is the lightest boson that can be exchanged by incoming 
neutralinos. Our corrections should then remain perturbative for WIMP masses up to at least 
1.2 TeV. 



4.3 Results and Discussion 

First we consider the scenario with wino-like LSP. This can e.g. be motivated from scenarios 
with anomaly mediated supersymmetry breaking, where the gaugino masses are related by 

12 Over the range shown in Fig. [7] the cross section strictly speaking does not violate unitarity, i.e. the 
annihilation amplitude is still smaller than unity. However, the behavior of the cross section at large LSP mass 
is clearly pathological, and would indeed lead to true unitarity violation at sufficiently large mass. 

13 Imaginary parts were introduced to regularize the infrared (IR) divergence in reactions like xt + Xi ~ ► 
7 + H + , which occurs for s ~ m #+ • In this special situation the photon energy is very small, so the ex- 
changed chargino is nearly on-shell. Introducing an imaginary part for the propagator of this nearly on-shell 
chargino indeed regularizes this divergence; however, the proper treatment of IR divergences instead requires 
the calculation of IR divergent one-loop diagrams, leading to an IR finite total result. Since this IR problem is 
relevant only for very special parameter choices, we kept the original regularization of DarkSUSY for final states 
containing one scalar and one massless gauge boson. 
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Figure 7: Annihilation rate va(xtXi ~>* for three-momentum \p\ = 10 _3 m^o. The 

solid (black) curve shows the prediction of the original DarkSUSY (version from November 2012), 
while the dot-dashed (red) curve, which shows the expected mZg behavior at large masses, is 
the prediction of the corrected version of DarkSUSY; see the text for details. 



M 2 ~ ijlfi . (53) 

We consider a wino mass between 100 GeV and 1.4 TeV. The higgsino and sfermion masses 
are set very high (30 TeV and 6 TeV, respectively), so that higgsino and sfermion exchange 
diagrams are very strongly suppressed. We assume that there is no flavor mixing in the sfermion 
sector. Not surprisingly, the spectrum calculator of DarkSUSY gives a light neutralino mass 
range from 100.0 GeV to 1.4 TeV. The original DarkSUSY code only calculates the spectrum 
up to the tree level, which underestimates the chargino-neutralino mass splitting in this case. 
Because the Sommerfeld correction is sensitive to the mass splitting, we add 0.17 GeV by hand 
to the chargino masses in the code. 

In this scenario, the (co-) annihilation processes that are relevant for the calculation of the 
LSP relic density are 



X° + xl^X + Y, 

Xi + Xi — ► X + Y, and its C-conjugate 
Xi + Xi — * X + Y, and its C-conjugate 

where X and Y stand for generic standard model particles. Sommerfeld-enhanced W exchange 
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can mix the first and second types of processes, whereas reactions of the third and fourth types 
only receive diagonal corrections, since the total charge in the initial and intermediate state 
must be the same. 

Here we focus on final states with sizeable annihilation rates a^v and plot the ratios of the 
corrections due to various intermediate states and the tree-level annihilation rates, 

R^-, (54) 
a 

against the WIMP mass. These corrections can be suppressed because of either of the following 
two reasons. It could be that a given intermediate state is only accessible via suppressed 
fermion-fermion-boson couplings. Examples are all process Xi + Xi ~~ Xi + Xi ~^ X + Y, 
where the rescattering Xi + Xi ~ >* Xi + Xi can De mediated by the Z boson or one of the 
CP-even neutral Higgs bosons. However, for a pure wino LSP the XiXi^ coupling is absent for 
the same reason that the SM doesn't have a triple— Z coupling, and the XiXi(h,H) couplings 
are absent because they require a non- vanishing higgsino component of Xv Numerically we 
find R < 10~ 6 in these cases. 

Small corrections also result if the given intermediate state has a very small annihilation 
cross section into the final state under consideration. For example, the correction to Xi + Xi ~~ ^ 
u + u annihilation from the x? + Xi intermediate state is very small since x? + Xi ~> u + u is 



suppressed by the very large m s mass we are considering: recall that XiXi has very suppressed 



couplings to both the Z and the neutral Higgs bosons, so that the s-channel contributions to 
Xi + Xi — ► u + w are also ver Y small. In contrast, xt + Xi u + u has sizeable tree-level 
annihilation rate since 7 and Z exchange in the s-channel contribute with full gauge strength. 
Again we find R < 10~ 6 in this case. Corrections of this size are obviously negligible. 

If neither of these two conditions is satisfied, corrections become quite large for large m^o 
and small three-momentum \p\ . As examples we show in Fig. [8] corrections to annihilation into 
W pair final states for different combinations of initial and intermediate states, and for two 
values of the three-momentum in the initial state in units of the LSP mass. Note that these 
figures only include corrections due to the exchange of massive bosons, in particular W ± and 
Z exchange. The XiXi (XiXi") initial states also contain (diagonal) Sommerfeld corrections 
due to photon exchange, which have been computed some time ago [21]. The corresponding 
exact (all-order) corrections are incorporated in our code, but have been suppressed "by hand" 
when producing the results shown in Fig. [8] in order to show more clearly the effect of one-loop 
corrections mediated by massive bosons. These photonic all-order corrections will be included 
later for the calculation of the relic density. 

Some features of these plots need to be remarked upon. First, in Fig. [8^l the three- 
momentum is always smaller than the mass of the exchanged boson (0 = W or Z). As a 
result, the correction is 0(am ? o/mA and hence increases monotonically with increasing LSP 
mass. On the other hand, in Fig. pb the three-momentum can become bigger than m^. The 
scale of the correction is then set by am%o/\p\, which is independent of the LSP mass in Fig. [8] 
since \p\ is taken to be a fixed fraction of the LSP mass here. In this case the corrections satu- 
rate beyond some value of the LSP mass. This is in accordance with the discussion of Section 

El 



14 Even for smaller the S-wave contribution to this cross section would be suppressed by a factor 
(m u /m^o) 2 . 



24 



1.5 - 



1 - 



Pi 0.5 - 



- 



-0.5 - 




-1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 

100 200 300 400 500 600 700 800 900 1000 1100 1200 1300 1400 1500 

WIMP mass [GeV] 



(a) \p\ — O.Olm^o 




-0.15 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 

100 200 300 400 500 600 700 800 900 1000 1100 1200 1300 1400 1500 

WIMP mass [GeV] 



(b) \p\ = 0.33m^o 

Figure 8: Relative size of corrections to annihilation into W pairs from different combinations 
of wino-like initial and intermediate states: X1X1 ~> X1X1 ~^ W + W~ (solid, black); XiXi ~^ 
X1X1 ~^ W + W~ (dot-dashed, red); and XiXi ~^ XiXi ~^ W + W + (dashed, blue). The upper 
(lower) frame is for cms three-momentum \p\ = O.Olm^o (\p\ = 0.33777,^0). 
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Secondly, while both annihilation reactions into W + W~ pairs receive positive corrections, 
the cross section for the annihilation of two positive charginos gets a negative correction because 
the potential between them is repulsive. Note that initial states containing two identical Dirac 
fermions, rather than a fermion-antifermion pair, lead to "clashing Dirac arrows" in our basic 
diagram of Fig. [2j As discussed in Sec. 2.1, we treat this using Denner's convention [TT], which 
gives an explicit minus sign in front of our basic correction of Eq. ^ . 

We finally note that for large LSP mass the relative correction to cf{xiXi ~~ W + W~) 
from the XiXi intermediate state is about two times bigger than the relative correction to 
a (XiXi ~^ W + W~) from the X1X1 intermediate state. The Xi ~ Xi mass splitting 5m is 
negligible even for |p| = O.Olm^o, because m^oSm <C m^,; see Eq.(34). Since both rescatterings 
proceed via W ± exchange, the Sommerfeld correction factors describing the rescattering are 
nearly the same in both cases. Moreover, while the XiXi initial state can form a spin triplet 
(total spin S = 1), the XiXi state has to be a spin singlet (S = 0) in the S'-wave. This means that 
the XiXi intermediate state can only give an S'-wave correction to the spin-singlet component 
of the XiXi annihilation reaction. In combination, these two facts imply that the numerators 
of R(xiXi -> XiXi -> W + W~) and R{xtxi ~> XiXi W + W~) are essentially the same. 
However, the correction to XiXi annihilation receives an extra factor of 2 (in the 5- wave), since 
the rescattering can produce both a XiXi an d XiXi intermediate state, whereas the XiXi 
intermediate state is unique; see the discussion in Sec. 2.1. On the other hand, the denominators 
of the two corrections R are quite different, since cr(XiXi ~~ >" W + W~) ~ ^{XiXi ~^ W + W~). 
This over-compensates the relative factor of 2 between the two numerators, making the relative 
correction to chargino annihilation about two times larger than that for neutralino annihilation. 
Moreover, chargino annihilation also receives sizable "diagonal" corrections from the chargino 
pair intermediate state, accessible via Z and 7 exchange; since no boson has sizable diagonal 
couplings to two neutralinos in this case, the "diagonal" corrections to LSP annihilation are 
negligible in this example. 

In the next step, we use the corrected annihilation rates to calculate the relic density for 



LSP mass between 100 GeV and 1.5 TeV. The result is plotted in Fig. [9aJ which also shows the 
tree-level prediction. Evidently the increase of the (co-) annihilation cross sections reduces the 
relic density. Note that within the range of x masses where the correction remains perturbative 
even for small three-momentum, the thermal relic density in standard cosmology (which is 
assumed here) is well below the total required Dark Matter density. This problem can be 
solved by introducing a second Dark Matter component, e.g. an axion or axino [25]. Another 
possibility is to enhance the expansion rate of the Universe (i.e., the Hubble parameter) during 
Xi decoupling, which increases the thermal Xi renc density [26J; in that case the relative size 



of the correction to the relic density would still be similar to that shown in Fig. 9a 



This relative correction, SQ^o/Q^o, is plotted in Fig. 9b We see that the one-loop Sommer- 
feld correction can be as high as 20%. As expected, the corrections become more important at 
higher LSP mass. The curve flattens towards large LSP mass, since thermal averaging of the 
co-annihilation cross sections favors three-momenta \p\ ~ 0.1 to 0.2m^,o, which becomes larger 
than the mass of the exchanged boson = My/ or Mz for large LSP mass; this is the same 
effect we saw (for slightly higher three-momentum) in Fig. |8j 

Next we turn to scenarios with higgsino-like LSP. For simplicity we assume that gaugino 
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Figure 9: Thermal LSP relic density with and without the one-loop "Sommerfeld" correction 
for wino-like LSP. The right panel shows the relative size of the correction to the relic density. 
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masses unify, which implies for weak-scale masses: 

Mx = ^M 2 . (55) 

In practice this does not matter, since all gaugino masses are set very high, M\ = 9.5 TeV. We 
consider higgsino masses between 100 GeV and 1.4 TeV, which leads to a very similar range for 
the LSP mass. Due to the very large gaugino masses, the (tree-level) mass splittings between 
the higgsino-like states amounts to at most 0.6 GeV. The sfermions are again assumed to be 
very heav}|^| (15 TeV), and the flavor mixing in the sfermion sector is turned off. 

As mentioned in the beginning of this Section, this situation is somewhat more complicated 
to analyze than scenarios with wino-like LSP, since there are three distinct higgsino-like states, 
but only two wino-like ones. Correspondingly, the following initial and intermediate states have 
to be considered: x°iX°i, XiXf, X^xf, X1X1] XiXv X1X2, XiXi- Since higgsino-gaugino mixing 
is again very small, Higgs boson exchange corrections can again be neglected. However, Z 
boson exchange corrections are now also sizable for neutralino initial states, since the Z\\\\ 
coupling is large, although the diagonal Zy^xl, {i — 1)2) couplings remain small. 

As before, a certain intermediate state only leads to a significant correction to a given 
annihilation reaction if both the rescattering rate from the initial to the intermediate state and 
the annihilation rate from the intermediate to the final state are large. Moreover, we again find 
that for fixed velocity of the annihilating particles, the size of the corrections increases with 
the LSP mass unless the three-momentum in the initial state is much larger than the mass of 
the exchanged boson. 

One important, qualitatively new feature emerges in this scenario: negative corrections 
become common, not restricted to annihilation processes of fermions with the same charge. 
This is possible because of the relative phase between the lightest and the next-to-lightest 
neutralino. 

To see this, consider the limit where Xi an d X2 are pure higgsino states. For positive \i the 
neutralino mixing matrix Z now has the forrrj^l 

/ l/y/2 -l/V2\ 
z= i/V2 i/V2 (56) 

\ / 

In addition, the chargino mixing matrices U and V take the form 

K-v-gJ). (57) 

From these we can calculate the relevant couplings between \ states and electroweak gauge 
bosons, using the Feynman rules listed in the Appendix of ref. [TU]; for the convenience of the 
reader we include the relevant rules in Appendix A. It is then not difficult to see that negative 
corrections can appear in many cases. 



15 This means that fermion-sfcrmion loop contributions to the higgsino mass splitting |14U15) . which we have 
ignored, will also be small. 

16 For /i < the symmetric state is the lighter one, and the i still appears for the second neutralino; i.e. 
Z13 = Z14 = —iZ23 = iZ2i = l/v2- Moreover a sign then appears in one of the non-vanishing entries of either 
U or V. The subsequent discussion still goes through in this case. 
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Figure 10: An example where the total correction can be negative. Note that in the right 
diagram, one of the intermediate Xi has been replaced by this gives a relative minus sign 
between the contributions of these two diagrams. 

As a first example, consider corrections to xt + Xi ~~ W + + W~ annihilation via neutralino 
pair intermediate states. Specifically, compare the XiX? an d X1X2 intermediate states with 



t-channel annihilation, as shown in Fig. 10 (The X1X1 P & i r basically does not annihilate via 
an s-channel diagram.) The relative sign between these contributions depends on the products 
of couplings associated with the four vertices shown in Fig. 10; in contrast to Eq.Q, Tj now 



includes both the coupling strength and the Dirac structure of the corresponding vertex. 

Let us begin with the X1X1 intermediate state. The two vertices F±, T 2 describing the 
rescattering stage are: 



Tt : ig^{C[ x P h + C*P R ) = ig 2 Y 
T 2 : ig 2l v (C^P L + C$P R ) = ig 2l v 



Similarly, the two vertices T 3 , T 4 describing the t-channel annihilation of the intermediate state 
are: 

1 

2 



T 3 : igzYiC^PL + C? X P R ) = ig 2l p ■ 
T 4 : i92-f{C%P L + C^P R ) = ig 2l c 
Now consider the X1X2 intermediate state. The four corresponding vertices are: 

0' 

1 

1— 
2 

I) 

1 

1— 
2 





igz-fV&PL + C*P R ) -- 


= i027" • 




tg 2 Y(C 2 L 1 *P L + C**P R ) 


= W2Y 


r' 
1 3 


tg 2 Y(C^P L + CKp R ) = 


= W2Y ■ 


r' 
1 4 


ig 2 -f(C£P L + C*P R ) 


= W2l U 
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From these expressions we see that the two intermediate states lead to exactly the same Dirac 
structure of the vertices [all are vector-like, because the two higgsino doublets form a vector-like 
representation of SU(2)]; moreover, all couplings have the same strength (i.e., absolute value). 
However, due to the two i factors appearing for the X1X2 intermediate state, the product of 
couplings for this intermediate state is negative, while it is positive for the intermediate 
state. 

A similar analysis shows that the combination of coupling factors for the X2X2 intermediate 
state is the same as that for the X1X1 intermediate state, while the X2X1 intermediate state 
contributes with the same product of couplings as the X1X2 intermediate state. In the limit 
where all mass differences between the charged and neutral higgsinos can be ignored, the total 
contribution from all four intermediate states thus vanishes! 

This may be surprising at first sight, but it can be understood from the observation that 
in this limit, the higgsinos form a degenerate SU(2) doublet of Dirac fermions. In that case a 



diagram like those in Fig. 10 does not exist. Start with the incoming xt ■ It can emit a W + 
at the first vertex to turn into the neutral Dirac higgsino, which is the lower component of the 
doublet. However, this lower component cannot emit yet another W + at the second vertex, as 
required in these diagrams, so they do not exist in this limit. 

On the other hand, the neutral Dirac higgsino can emit a W~ , so a w-channel diagram similar 



to the ones shown in Fig. 10 should exist in the pure higgsino limit. In our calculations with two 
distinct neutral Majorana higgsinos, the coupling factors for these ^-channel diagrams can be 
obtained from those for the t-channel by swapping the couplings C in T% and T^, taking care to 
keep track of neutralino indices. This does not change anything for the X1X1 intermediate state, 
where the same neutralino index appears everywhere. However, for the X1X2 intermediate state, 
we now get the coupling factors Cji, in T' 3 without complex conjugation; instead, now C[* 
and appear in T^, but these couplings are real. The purely imaginary couplings therefore 
now contribute i x (— i) = 1, rather than i x i — — 1 for the t-channel diagrams. Hence there 
is no cancellation between the different w-channel diagrams] 17 ] 



Fig. 11 shows that the total correction to X1X1 annihilation from neutralino intermediate 
states is nevertheless negative. The reason is that only intermediate states containing two 
different neutralinos can annihilate efficiently into SM fermion-antifermion final states via s- 
channel exchange of a Z boson. The total contribution from these mixed intermediate states is 
therefore considerably larger in magnitude than the contribution from both intermediate states 
containing two equal neutralinos. The fact that the X1X2 intermediate state gives a negative 
correction to cr(xiXi ~~ >* ff)i where / is an SM fermion, can again be understood easily using 
the notion of Dirac higgsinos: the coupling of the Z boson to neutral and charged higgsinos will 
then have opposite sign, since their J 3 (weak isospin) values have opposite sign. Note, however, 
that there are also positive contributions from X1X1 intermediate states, which are not shown 



in Fig. 11 



Another instance that involves a possibly negative correction occurs in X1X1 annihilation. 
We compare the intermediate states X1X1 an d X1X21 see Fig- 12 We only consider s-channel 



annihilation in this example. The left diagram then shows the only sizable contribution from 
X1X1 ( or y i ce versa) intermediate states. We again calculate the three vertex factors in these 



17 The Dirac higgsino limit also allows to understand why a(xtxt ~^ W + W + ) is very small for higgsino-like 
LSP, in sharp contrast to the case of wino-like LSP where it is large. 



30 



0.05 



Pi 



-0.05 



-0.1 



-0.15 



-0.2 



intermediate x"i>X U i = 
intermediate X°1.X°2 



100 200 300 400 500 600 700 800 900 1000 1100 1200 1300 1400 1500 

WIMP mass [GeV] 



Figure 11: The relative correction to the total annihilation rate of XiXi-, f° r the X1X1 (solid, 
black) and X1X2 (dashed, red) intermediate states; the contribution from X2X2 i s very similar 
to that from XiXv The contributions from all Standard Model final states have been summed, 
and \p\ = O.Olm^o. 



diagrams. For the left diagram: 

Ti = ig 2 ^(C^P L + C*P R ) = igrf 
T 2 = ig 2 r{C?*P L + C£P R ) = tg 2 Y 
To = igzlUC&PL + C^P R ) = igrf 



For the right diagram: 



r; = %g % rf{c\ x p L + c*p r ) = i&rf 

T' 2 = igxf{C%P L + C^Pr) = ig 2 Y 



= ig*y»{C*P L + C^Pr) = ig 2 <f ■ 



2) 



From Eq.(52) one again finds a relative minus sign between these contributions. In fact, in the 



pure higgsino limit these contributions will cancel exactly; this can also be understood from 
the observation that no such diagram can be drawn for neutral Dirac higgsinos. 
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Figure 12: A second example where the total correction can be negative. In the right dia- 
gram the intermediate Xi nas been replaced by a giving a relative minus sign between the 
contributions from these diagrams. 




Figure 13: The relative correction to the total annihilation rate of xTxi> for the XiX.t (solid, 
black) and X2X1 (dashed, red) intermediate states. The contributions from all Standard Model 
final states have been summed, and \p\ = O.Olm^o. 



We see in Fig. 13 that after summing over all SM final states, the negative contributions 
again win. In this case the two intermediate states shown have very similar annihilation cross 
sections for all contributing final states. However, the intermediate state containing the heavier 
neutralino x\ * s enhanced because it is also accessible via Z exchange in the rescattering process, 
whereas the X1X1 intermediate state is only accessible via the W exchange diagram shown in 



Fig. 12 



In fact, the numerical calculation shows that negative corrections occur quite frequently in 
the co-annihilation of higgsino-like states. For example, X1X1 states coupled to X1X2 states, 



32 



with either state in the initial and the other in the intermediate state, yields negative corrections. 
The same is true for XiXi coupling to X1X21 where again either state can be in the initial state. 

The cancellation of positive and negative corrections to individual annihilation processes in 
the end leads to a small correction to the total effective annihilation rate, which can be either 
positive or negative. As a result, the relic density is slightly enhanced or reduced over the range 
of the WIMP mass considered, as shown in Fig. 14b, For the correct relic density (Qh 2 ~ 0.113), 
the total one-loop correction turns out to be very small, less than 0.5%. However, this is largely 
"accidental", since the corrections for fixed initial, intermediate and/or final state are typically 
much larger, as shown in the previous figures. 



5 Summary and Conclusions 

In this paper we have computed enhanced one-loop corrections to the co-annihilation of WIMPs 
due to the exchange of a light boson in the initial state, treating both S- and P-wave initial 
states and carefully including the effects from multiple interfering intermediate states (so-called 
"Sommerfeld" corrections). The ultimate goal is a more accurate calculation of the thermal 
WIMP relic density. 

In Sec. 2 we extended the formalism of ref. [5] to deal with the "multistate" Sommerfeld 
effect, where the particles in the intermediate state could be different from those in the initial 
state, although the mass splitting should be relatively small. Since co-annihilation of the 
slightly heavier partners of the WIMPs also needs to be treated, we considered cases where the 
intermediate state is lighter or heavier than the initial state, in addition to the usual case where 
the initial and intermediate states have the same masses. We found exact analytical expressions 
for the functions describing the one-loop corrections for all three cases; these supercede the 
numerical fits found in ref. [3] for the case of equal masses. 

As the intermediate state particles are almost on-shell, the boson exchange can still be 
regarded as a rescattering reaction, which however is in general off-diagonal. As a result, the 
corrections no longer factorize at the level of the annihilation cross sections, although they do 
factorize at the level of the amplitude. The existence of several interfering intermediate states 
can lead to additional complications, as discussed in Sec. 2.1. In the final subsection of Sec. 2 
we showed that the exchange of a light fermion does not lead to enhanced corrections to the 
co-annihilation of a boson with a (heavy) fermion. 

The dependence of the loop functions on various quantities is discussed in Sec. 3. We 
found that the mass splitting 5m x between co-annihilating particles affects the loop functions 
significantly whenever m x \5m x \ > m"^, where m x is the WIMP mass and m$ is the mass of 
the exchanged boson. For very small external three-momentum a non-vanishing mass splitting 
always reduces the correction, the effect being more pronounced for anihilation from the S- 
wave. However, if the intermediate state is heavier than the initial state, the loop function 
develops a peak where the center-of-mass frame energy equals exactly the total mass of the 
intermediate state. 

In Sec. 4 we applied this formalism to the calculation of the relic density of the lightest 
neutralino in the MSSM. In that case co-annihilation is generic if the LSP is either wino- 
or higgsino-like. In the former case the co-annihilation with the lightest chargino has to be 
considered. We found that most corrections are positive, i.e. they reduce the relic density 
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(b) relative correction to Vlh 2 

Figure 14: Current relic density with and without the one-loop "Sommerfeld" correction for 
higgsino-like LSP. The right panel shows the relative size of the correction. 



even further. The correct thermal relic density is then reached for a range of WIMP masses 
where the one-loop corrections become so large that they need to be re-summed [2]. For 
higgsino-like LSP, the two lightest neutralinos and the lighter chargino all contribute in various 
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co-annihilation reactions. In this case many corrections turn out to be negative. We saw that 
in many cases this can be understood in the limit of exact higgsino LSP, in which case the two 
lightest neutralinos can be grouped into a neutral Dirac higgsino, which is an SU(2) partner 
of the lighter chargino. We found that in this case the total correction to the thermal relic 
density happens to cancel to good approximation, even though corrections to specific initial 
and/or final states can be quite sizable. 

This paper thus completes the model-independent treatment of one-loop "Sommerfeld"- 
enhanced corrections to WIMP annihilation, and at the same time adds to the growing literature 
on potentially large corrections to the (co-) annihilation of supersymmetric neutralinos. 
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Here we list some Feynman rules for vertex factors in the MSSM used in our calculations. They 
are taken from Ref.[10j and adapted to Denner's convention 



A MSSM Vertices 
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Figure 15: x + X°W + vertices. 
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Figure 16: X + X + l vertex. 
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Figure 17: vertex. 
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Figure 18: x x Z vertex. 
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